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For compact U(l) lattice gauge theory (LGT) we have performed a finite size scaling analysis on 
N T Ns lattices for N T fixed by extrapolating spatial volumes of size N s < 18 to N s — » oo. Within 
the numerical accuracy of the thus obtained fits we find for iV T = 4, 5 and 6 second order critical 
exponents, which exhibit no obvious N T dependence. The exponents are consistent with 3d Gaussian 
values, but not with either first order transitions or the universality class of the 3d XY model. As 
the 3d Gaussian fixed point is known to be unstable, the scenario of a yet unidentified non-trivial 
fixed point close to the 3d Gaussian emerges as one of the possible explanations. 
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I. INTRODUCTION 

Abelian, compact U(l) gauge theory has played a 
prominent role in our understanding of the permanent 
confinement of quarks. It was first investigated by Wil- 
son in his 1974 milestone paper [1], which introduced 
lattice gauge theory (LGT). For a 4d hypercubic lattice 
his U(l) action reads 

S({U}) = J2Sn (1) 

□ 

with S n = Re (Ui^Uj^Ui^Uj^), where i\, ji, i 2 and 
j2 label the sites circulating about the square □ and 
the Uij are complex numbers on the unit circle, £/y = 
exp(i4>ij), < fa < 2ir. 

Wilson concluded that at strong couplings the the- 
ory confines static test charges due to an area law for 
the path ordered exponentials of the gauge field around 
closed paths (Wilson loops). A hypothetical mechanism 
of confinement was identified by Polyakov [2], who at- 
tributed it in 3d Abelian gauge theory to the presence of 
a monopole plasma. For the 4d theory at weak coupling 
both Wilson and Polyakov expected a Coulomb phase in 
which the test charges are not confined. The existence of 
two distinct phases was later rigorously proven [3] . 

So it comes as no surprise that 4d U(l) LGT was the 
subject of one of the very early Monte Carlo (MC) cal- 
culations in LGT [4]. One simulates a 4d statistical me- 
chanics with Boltzmann factor exp [—/3 g S({U})] and pe- 
riodic boundary conditions (other boundary conditions 
are possible too, but are not considered here), f3 g = 1/g 2 
is related to the gauge coupling g 2 , (3 g — is the strong 
and /3 g —> oo the weak coupling limit. The study [4] al- 
lowed to identify the confined and deconfined phases. Af- 
ter some debate about the order of the phase transition, 
the bulk transition on symmetric lattices was suggested 
to be (weakly) first order [5] , a result which was substan- 
tiated by simulations of the Wuppertal group [6,7] . Other 
investigations followed up on the topological properties 
of the theory. This lies outside the scope of the present 



paper. The interested reader may trace this literature 
from [8]. 

The particle excitations of 4d U(l) LGT are called 
gauge balls and in the confined phase also glueballs. 
Their masses were first studied in Rcf. [9]. In the con- 
fined phase all masses decrease when one approaches the 
transition point. Crossing it, they rise in the Coulomb 
phase with exception of the axial vector mass, which is 
consistent with the presence of a massless photon in that 
phase. Recently this picture was confirmed in Ref. [10], 
relying on far more powerful computers and efficient noise 
reduction techniques [11]. The first order nature of the 
transition prevents one from reaching a continuum limit, 
as is seen in Fig. 7 of [10]. In contrast to that investiga- 
tions in a spherical geometry [12] and of an extended U(l) 
Wilson action [13] reported a scaling behavior of glueballs 
consistent with a second order phase transition. But this 
is challenged in other papers [14,15], so that it remains 
questionable whether an underlying non-trivial quantum 
field theory of the confined phase can be defined in this 
way. 

Here we focus on U(l) LGT in finite temperature ge- 
ometries. We consider the Wilson action (1), choose units 
a = 1 for the lattice spacing and perform MC simula- 
tions on N T Ns lattices. Testing U(l) code for our biased 
Metropolis- heatbath updating (BMHA) [16], we noted 
on small lattices that the characteristics of the first or- 
der phase transition disappeared when we went from the 
N T = N s to a N T N^, N T < N s geometry. This moti- 
vated us to embark on a finite size scaling (FSS) calcula- 
tion of the critical exponents of U(l) LGT in the N T N%, 
N T = constant, N s — > oo geometry. For a review of FSS 
methods and scaling relations see [17]. 

Later we learned about a paper by Vettorazzo and de 
Forcrand [18], who speculate about a scenario of two 
transitions at finite, fixed N T : One for confinement - 
deconfinement, another one into the Coulomb phase, 
both coinciding only for the zero temperature transition. 
Their claim for the confinement-deconfmement transition 
is that it is first order for N T = 8 and 6, N s — > oo, be- 
coming so weak for N T < 4 that it might then be second 
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order. In contrast to having two transitions at finite N T 
the conventional expectation appears to us one transi- 
tion, which is second order and in the 3d XY universality 
class, switching to first order for sufficiently large N T . See 
Svetitsky and Yaffe [21] for an early discussion of some 
of these points. 

In the next section we present our numerical results 
in comparison with previous literature, followed by sum- 
mary and conclusions in the final section. 

II. NUMERICAL RESULTS 

Our FSS analysis relies on multicanonical simulations 
[19] for which the parameters were determined using a 
modification of the Wang-Landau (WL) recursion [20]. 
A speed up by a factor of about three was achieved 
by implementing the biased Metropolis-Heatbath algo- 
rithm [16] for the updating instead of relying on the 
usual Metropolis procedure. This is substantial as, for 
instance, our 16 4 lattice run takes about 80 days on a 
2 GHz PC. Additional overrelaxation [22] sweeps were 
used for some of the simulations. 

Our temporal lattice extensions are N T = 4, 5 and 6. 
For N s our values are 4, 5, 6, 8, 10, 12, 14, 16 and 18. 
Besides we have simulated symmetric lattices up to size 
16 4 . The statistics analyzed in this paper is shown in 
table I. The lattice sizes are collected in the first and 
second column. The third column contains the number 
of sweeps spent on the WL recursion for the multicanon- 
ical parameters. Typically the parameters are frozen af- 
ter reaching / = e 1 / 20 for the multiplicative WL factor 
(technical details of our procedure will be published else- 
where). Column four lists our production statistics from 
simulations with fixed multicanonical weights. Columns 
five and six give the (3 values between which our Markov 
process cycled. Adapting the definition of chapter 5.1 of 
[23] one cycle takes the process from the configuration 
space region at /3 m i n to (3 max and back. Each run was 
repeated once more, where after the first run the multi- 
canonical parameters were estimated from the statistics 
of this run. Columns seven and eight give the number of 
cycling events recorded during runs 1 and 2. 

Using the logarithmic coding of chapter 5.1.5 of [23] 
physical observables are reweighted to canonical ensem- 
bles. Error bars as shown in figures are calculated using 
jackknife bins (e.g., chapter 2.7 of [23]) with their num- 
ber given by the first value in column four (always 32), 
while the second value was also used for the number of 
equilibrium sweeps (without measurements) performed 
after the recursion. Weighted by the number of their 
completed cycles, the results from two or more runs are 
combined for for the final analysis (compare chapter 2.1.2 
of [23]). 



TABLE I. Statistics of our MC calculations. The simula- 
tion with * attached in the WL column uses 22 WL recursions, 
all others 20. 
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A. Action variables 

Figures 1 and 2 show for various values of N s the spe- 
cific heat 

= ^[(S 2 )-(S) 2 } with N = N T N! (2) 

in the neighborhood of the phase transition for N T = 6 
and on symmetric lattices. The (3 ranges in the figures 
are chosen to match. 

In Fig. 3 we show all our specific heat maxima on a 
log-log scale. Our data for the symmetric lattices are for 
N s > 8 consistently described by a fit to the first order 
transition form [24] C max (N s ) / (6N) = c G +a 1 /N+a 2 /N 2 . 
The goodness of our fit is Q = 0.64 (see, e.g., chap- 
ter 2.8 of Rcf. [23] for the definition and a discussion 
of Q), and its estimate for the specific heat density is 
c = 0.0001961 (26). This is 10% higher than the c value 
reported by the Wuppertal group [7], where lattices up 
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FIG. 1. Finite size dependence of the specific heat func- 
tions C(/3) on N T — 6 lattices. 
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Maxima of the specific heat. 
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FIG. 2. Finite size dependence of the specific heat func- 
tions C(/3) on iV T = N s lattices. 

to size 18 3 were used. For the interface tension consistent 
fits to a small positive as well as to a zero infinite volume 
value are possible. Interestingly the Wuppertal group de- 
cided against including their largest lattices from [6] into 
their final analysis [7], because action data on them did 
not cover the double peak region coherently (communi- 
cated by T. Neuhaus). 

For N T — 4, 5 and 6 our curves in Fig. 3 are linear fits 
in N s , C max (N s ) = ai N s +a + a_i/iV s . For N T = 4 the 
goodness of this fit is Q = 0.20 using our N s > 6 data. 
But for N T — 5 and 6 the Q values are unacceptably 
small, although the data scatter nicely about the curves. 
For large N s the maxima of the specific heat curves scale 
like (see [17]) 



Cmax(^s 



NT /v ■ 



(3) 



terns become three-dimensional, so that a/v = 3 would 
be indicative of a first order transition, while our data 
are consistent with the second order exponent a/v = 1. 

This has to be contrasted with the claim by Vettorazzo 
and de Forcrand [18] that the N T > 6 transitions are first 
order. For iV r = 8 and 6 their evidence relies on simu- 
lations of very large lattices. Differences in action values 
obtained after ordered and disordered starts support a 
non-zero latent heat in the infinite volume limit. For 
N T = 6 the spatial lattice sizes used are N s — 48 and 60 
and their MC statistics shown consists of 5 000 measure- 
ments per run, separated by one hcatbath plus four over- 
relaxation sweeps (these units are not defined in [18], 
but were communicated to us by de Forcrand and previ- 
ously used in [15]). For a second order phase transition 
the integrated autocorrelation time T; n t scales approxi- 
mately ~ N% and we estimate from our own simulations 
on smaller lattices that in units of those measurements 
Tint « 7 000 for N T = 6 and N s = 48. A MC segment 
of the length of Tint delivers one statistically independent 
event (e.g., chapter 4.1.1 of [23]). Therefore, the run of 
[18] would in case of a second order transition be based 
on less than one event and strong metastabilities would 
be expected as soon as the Markov chain approaches the 
scaling region. For N s = 60 and the N r — 8 lattices 
the situation is even worse. We conclude that these data 
cannot decide the order of the transition. 

Let us remind the reader that a double peak alone 
does not signal a first order transition. One has to study 
its FSS behavior, but no error bars can be estimated 
when one has only one statistically independent event. 
Actually for our larger spatial volumes we find double 
peaks in our 6 x Nf action histograms and they are also 
well-known to occur for the magnetization of the 3d Ising 
model at its critical point [25]. 



where one has a/v = 4 in case of the first order transition 
for N T = N s . In the N T fixed, N s — > oo geometry the sys- 
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B. Polyakov loop variables 

Besides the action we measured Polyakov loops and 
their low-momentum structure factors. For U(l) LGT 
Polyakov loops are the Uij products along the straight 
lines in N T direction. Each Polyakov loop Pg is a com- 
plex number on the unit circle, which depends only on 
the space coordinates, quite like a XY spin in 3d. We 
calculate the sum over all Polyakov loops on the lattice 



(4) 



The critical exponent 7/1/ is obtained from the maxima 
of the susceptibility of the absolute value \P\, 
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Structure factors are defined by (see, e.g., Ref. [26]) 
^ P(f) exp(ifcr) 



(5) 
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where n is an integer vector, which is for our measure- 
ments restricted to (0,0,1), (0,1,0), and (1,0,0). Max- 
ima of structure factors scale like 



(8) 



The exponents can be estimated from two parame- 
ter fits (A) Y = a\ Ng 2 . Due to finite size corrections 
the goodness Q of these fits will be too small when all 
lattice sizes are included. The strategy is then not to 
overweight [27] the small lattices and to omit, starting 
with the smallest, lattices altogether until an acceptable 
Q > 0.05 has been reached. We found a rather slow con- 
vergence of the thus obtained estimates with increasing 
lattice size. This can improve by including more parame- 
ters in the fit. So we used the described strategy also for 
three parameter fits (B) Y = ao + ai Y" 2 . The penalty for 
including more parameters is in general increased insta- 
bility against fluctuations of the data and, in particular, 
their error bars. For a number of our data sets this is the 
case for fit B, so that an extension to more than three 
parameters makes no sense. We performed first the fit B 
for each data set, but did fall back to fit A when no con- 
sistency or stability was reached for a fit B including at 
least the five largest lattices. The thus obtained values 
are listed in table II. Table III gives additional informa- 
tion about the fits. 

Our lattices support second order transitions for Y T = 
4, 5 and 6. The evidence is best for observables derived 



TABLE II. Estimates of critical exponents as explained in 
the text. Properties of the fits are summarized in table III. 
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TABLE III. Number of data used and type of fit (A or B 
as explained in text), goodness of fit Q. 
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from Polyakov loops. For example, in Fig. 4 we show our 
data for the maxima of the Polyakov loops susceptibility 
together with their fits used in table II (for the sym- 
metric lattices the data are connected by straight lines). 
For fixed N T we find an approximately quadratic increase 
with N s , while there is a decrease for the symmetric lat- 
tices, which appears to converge towards zero or a finite 
discontinuity (note that one has no common scale for 
Polyakov loops from symmetric lattices, because their 
lengths change with N T ). 

Our structure factor data support that one is for 
(3 > p c in the Coulomb phase: As shown for N T = 6 in 
Fig. 5 the structure factors remain divergent for (3 > (3 C , 
as expected for a power law fall-off of Polyakov loop cor- 
relations. These observations apply to the (3 ranges (com- 
pare table I) covered by our multicanonical simulations. 
To have still reasonably many cycling events on large 
lattices, this range was chosen to shrink with increasing 
lattice size. So we test not very far into the /? > (3 C phase. 

The Polyakov loops describe 3d spin systems. So one 
would like to identify whether the observed transitions 
are in any of their known universality classes. At first 
thought the universality class of the 3d XY model comes 
to mind (e.g., [21]), because the symmetry is correct. It 
is easy to see that the N T = 1 gauge system decouples 
into a 3d XY model and a 3d U(l) gauge theory. The 
latter has no transition and is always confined. But one 
cannot learn much from this observation as there is no 
interaction between the two systems. Surprisingly the 
data of table II do not support the XY universality class. 
Although our estimates of 7/1/ agree with what is ex- 
pected, a/v is entirely off. For the XY model a small 
negative value is established [17], while Fig. 3 shows that 
all our specific heat maxima increase steadily. We re- 
mark that the scenario may change for N T < 4. We have 
preliminary results for N T — 2 and 3. The increase of the 
specific heat maxima becomes considerably weaker than 
for N T = 4. For N T = 2 it slows continuously down with 
increasing lattice size (so far up to N s = 20) and one can 
imagine that it comes altogether to a halt. Once com- 
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FIG. 4. Maxima of Polyakov loop susceptibilities. 



FIG. 6. Rescaled maxima of Polyakov loops susceptibili- 
ties. 
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FIG. 5. Structure factors (7) for N T = 6. 

pleted, our simulations for N T = 2 and 3 will be reported 
elsewhere. 



III. SUMMARY AND CONCLUSIONS 

In view of expected systematic errors due to our limited 
lattice sizes, one can state that our estimates of table II 
are consistent with the Gaussian values a/v = 1 and 
j/v = 2 (with error bars 0.3 for a/v and 0.1 for j/v). 
Using the hyperscaling relation 2 — a = dv with d = 3 
yields a = v = 1/2. The other estimates of exponents 
listed in table II provide consistency checks as they are 
linked to a/v = 1 and j/v = 2 by the scaling relations 
a + 2/3 + 7 = 2 and -f/v = 2 — 77. For the Gaussian expo- 
nents (1 — (3)/v = 1.5 and 77 = follows, both consistent 
with the data of the table. 

However, the problem with the Gaussian scenario is 
that the Gaussian renormalization group fixed point in 



3d has two relevant operators [28]. So one does not un- 
derstand why the effective spin system should care to 
converge into this fixed point [21]. Therefore, the inter- 
esting scenario of a new non-trivial (n-t) fixed point with 
exponents accidentally close to 3d Gaussian arises. An 
illustration, which is consistent with the data, is given in 
the last row of table II. The mean values are constructed 
to fulfill the scaling relations and match with v = 0.482, 
a = 0.554, 7 = 0.94, /3 = 0.253, 77 = 0.05. 

One may expect that the first order transition of the 
symmetric lattices prevails once N T is larger than the 
correlation length on symmetric lattices. But a non-zero 
interface tension has never been established for this tran- 
sition. So one could also imagine an instability under the 
change of the geometry. From a FSS point of view it ap- 
pears then natural that the character of the transition 
will not change anymore, once a value of -/V T has been 
reached, which is sufficiently large to be insensitive to lat- 
tice artifacts. Up to normalizations data from N T N^ and 
2N T (2N S ) 3 , N s > N T lattices should then become quite 
similar. We illustrate this here by rescaling the maxima 
of our Polyakov loop susceptibilities with a common fac- 
tor, so that they become equal to 1 on symmetric lattices. 
On a log-log scale the results are then plotted in Fig. 6 
against N S /N T . The behavior is consistent with assum- 
ing a common critical exponent for all of them (parallel 
lines are then expected for large N S /N T ). 

The litmus test for identifying a second order phase 
transition is that one is able to calculate its critical ex- 
ponents unambiguously. Instead of starting with data 
of uncontrolled quality from very large lattices, the FSS 
strategy is to control finite size effects by working the way 
up from small to large systems. With MC calculations 
FSS methods find their limitations through the lattice 
sizes, which fit into the computer and can be accurately 
simulated in a reasonable time. Within the multicanoni- 
cal approach "accurately" means that one has to get the 
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system cycling through the entire critical or first order 
region, and at least about one hundred cycles ought to 
be completed with measurements. 

Our lattice sizes are not small on the scale of typical 
numerical work on U(l) LGT, for instance the lattices 
used for the Wuppertal cq estimate of [7]. But we have 
not yet reached lattices large enough to provide hard ev- 
idence that there is no N s — > oo turn-around towards ei- 
ther a first order transition or the 3d XY fixed point. In 
particular in view of the fact that our data do not support 
the generally expected scenario, it would be desirable to 
extent the present analysis to the largest lattices that can 
be reached by extensive simulations on supercomputers, 
instead of relying on relatively small PC clusters. 

With mass spectrum methods [10,29] one may inves- 
tigate the scaling behavior of the model from a different 
angle. In particular observation of a massless photon 
[9] can provide more direct evidence for the Coulomb 
phase than our structure factor measurements. Finally, 
rcnormalization group theory could contribute to clarify- 
ing the issues raised by our data. Amazingly, even after 
more than thirty years since Wilson's paper [1] the na- 
ture of U(l) LGT phase transition is still not entirely 
understood. 
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